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Abstract 

An efficient Path Integral Monte Carlo procedure is proposed to simulate the behavior of quan- 
tum many-body dissipative systems described within the framework of the influence functional. 
Thermodynamic observables are obtained by Monte Carlo sampling of the partition function af- 
ter discretization and Fourier transformation in imaginary time of the dynamical variables. The 
method is tested extensively for model systems, using realistic dissipative kernels. Results are also 
compared with the predictions of a recently proposed semiclassical approximation, thus testing the 
reliability of the latter approach for weak quantum coupling. Our numerical method opens the 
possibility to quantitatively describe real quantum dissipative systems as, e.g., Josephson junction 
arrays. 
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The description of open quantum systems has significant imphcations in the fundamental 
concepts of quantum mechanics. It involves issues as the breakdown of the Schrodinger 
picture and the connection between the 'reversible' and the 'irreversible' world |Q. These 
basic problems gave the first strong boost to the research in this field in the early 50 's and 
still are among the most fascinating aspects of the topic. 

In the recent years a renewed interest in the study of quantum dissipation has come 
mainly from condensed-matter physics, as systems with an intermediate [mesoscopic) scale 
have been experimentally developed and theoretically analyzed. The fundamental issues in 
quantum and statistical mechanics are closely tied to technical applications, such as the 
single electron transistors 0, the Josephson quantum-bits [Q, the SQUIDS 0], and many 
others. In such systems, the characteristic quantum effects involve a macroscopic number 
of particles. The sizeable dimension of the devices implies that the relevant dynamical vari- 
ables can couple to a very large (infinite) number of degrees of freedom of the surrounding 
environment (or dissipation bath). The interaction between the open system and its envi- 
ronment leads in general to dissipation, fluctuation, decoherence and irreversible processes. 
These result in dramatic changes in the behavior of the system, for instance the dissipative 
phase transition in Josephson junction arrays 

At variance with classical thermodynamics, which is unaffected by dissipation, quantum 
thermodynamics is an ideal field to study the genuine interplay between quantum fluctua- 
tions and dissipation, which leads in general to interesting physics in the regimes of high 
quantum coupling and/or low temperature. Unfortunately, a suitable theoretical approach, 
allowing a faithful comparison with the experimental findings in the regime of high quantum 
coupling, is still lacking. In this Letter, we present an original Path Integral Monte Carlo 
(PIMC) approach, which allows one to efficiently deal with the quantum thermodynamics of 
many-body open systems in a wide range of couplings and temperature; we will show how 
working with Fourier transformed variables, and taking advantage of our knowledge of the 
exact quantum propagator of harmonic systems, allows us to get reliable results for many 
body-systems with reasonable numerical efforts. 

We address the problem of quantum dissipation within the Caldeira-Leggett(CL) frame- 
work [p, where dissipation is described as the result of a linear coupling of the physical 
system of interest with a bath of harmonic oscillators. By generalizing the CL formalism to 
many-body systems, the partition function of a quantum dissipative system is given by the 
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path integral 

Z = f V[q] e-'^'^K (1) 

with the Euchdean action 



S[q] 



du 



]^q{u)Aq{u) + V{q{u)) 



sAq] ■ (2) 



Here, q = {qi}i=i^...^M denotes the vector whose components are the M coordinates of the in- 
vestigated system and A = {Aij} is the mass matrix; the effects of dissipation are contained 
in the influence action, 

Sni[q] = 7^ / du du'q{u)K{u-u') q{u') . (3) 
Zii Jo Jo 

Within the many-body CL formalism the kernel K{u) = {Kij{u)} is an MxM matrix, 
which is in general nonlocal in space, i.e., nondiagonal. Thus the dissipation bath can drive 
also the spatial correlations of the system as it is expected, for instance, in the case of 
shunted Josephson junction arrays. The kernel matrix K{u), depends on the temperature 
T = is a symmetric and periodic function of the imaginary time u, K(u) = 

K{—u) = K{(3h — u), and has a vanishing average over a period. 

In order to numerically evaluate the integral appearing in Eq. (|l|), the standard PIMC 
method divides the imaginary-time interval [0,/?/i] into P slices of width e = f3h/P, P 
being the so called Trotter number. The coordinates q{u) turn into the discrete quantities 
Qi = 9(^^)5 P{ph)^^{qi — qi-i), and the partition function Z is obtained as the 

P — s> 00 extrapolation of 

Zp = CP-"^ n / dl^o / n dq^i e-^-l^''^^] , (4) 
1=1 •' i=\ 

where S'p[{q£}] is the discretized form of the action (g) and C is a temperature-independent 
normalization; the related macroscopic thermodynamic quantities are obtained through the 
accordingly generated estimators. 

The application of the standard PIMC approach to dissipative systems is made difficult 
by the fact that the kernel K{u—u') is usually not explicitly known in the imaginary-time 
domain, but rather in the Matsubara frequency one. In fact, it is generally given in terms 
of the bath spectral density or in terms of the Laplace transform of the damping function 
appearing in the phenomenological Langevin equation 0. This makes K{u—u') long-ranged, 
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cumbersome to evaluate and sometimes ill-defined as, e.g., in the case of the widely used 
Ohmic (or Markovian) dissipation. A possible way to avoid such problem could be the use 
of Fourier path-integral approaches, possibly supported by the partial-averaging scheme : 
the latter method, however, still deals with the problem of evaluating path averages of the 
potential over continuous paths. 

What we propose here is to start from the finite-P expression of the standard PIMC 
for the partition function and make there a lattice (discrete) Fourier transform, changing 
the integration variables from q^^ to qik by setting: 

qi = q+Y.Qke^'-''"' (5) 

A;=l 

SO that: 

M , , P-1 



M P-1 



exp <^ -J2 



k=l 



ph^ p" ' 2 



sin —A + -Kk 



Qk- 



~T.yi<ie)], (6) 



where Kk = {Kij^k} is the usual Matsubara transform Kk = du K{u) e'^"''^ of the 
dissipative kernel matrix at the Matsubara frequency Uk = 2'Kk/ (3h. 

We must however observe that the change of variables above is not enough to get a 
really efficient procedure to simulate many-body systems. In fact, in order to get a reliable 
thermodynamic limit finite-size effects have to be negligible; as a consequence the number M 
must be large enough, so that reaching high values of P may become computationally very 
demanding, making the extrapolation to P ^ cxd problematic. However, such difficulty can 
be largely circumvented by making use of our knowledge of both the finite- and infinite-P 
exact partition function of pure bilinear actions 0. According to Ref. [Eqs. (38) and 
(44)] any Monte Carlo estimate G{P) of a given quantity G obtained at finite P can be 
corrected by adding the exact (P — > oo) value Gha ^^^d subtracting the finite-P estimate 
GiS(P) of the same quantity |lO| for the (self-consistent) harmonic approximation of the 
dissipative action, getting: 

Gha(P) = G{P) + [Gffi - gS(P)] . (7) 

As it will clearly appear in the following applications, the last step turns out to be essential 
and truly effective in the investigation of many-body systems. 



As a preliminary test, we take a single particle of mass m in a non-linear potential 
V{q) = ef (g/cr), where e and a define the energy and length scale, respectively. In order to 
better examine the combined effects of quantum fluctuations and dissipation, it is useful to 
introduce the reduced temperature t = (/5e)~^ and the quantum coupling g = huo/e, where 
Uq = Jec^/ma^ and 



Xm being the absolute minimum of v{x). For odd Trotter 



number, P = 2N + 1, we get: 



p r r ^ f 1 

Zp = Ct^ dx Yl dakdbk exp <^ - — ^ v{xe) + 



N 

E 

fc=i 



— sm — + —Kk 

9^ P t . 



iai + 6 



(8) 



where x^ = x + 2X]fcLi (ofc cos + fe^sin^^), Kk = K^/uj^ and we have used the sym- 
metry properties of K{u), so that np-k = i^k- The real Fourier variables x, and are 
dimensionless, and the integrals in Eq. (||) may be numerically evaluated by standard Monte 
Carlo sampling techniques, e.g. the Metropolis one. 

Fig. shows the results obtained for the average potential energy when a quartic double- 



well potential v{x) 



X 



2\2 



and Ohmic dissipation, i.e., = 27i(t/ g)Tk, are considered 



(r is the damping strength in units of ujq)] the same model was already investigated in 
Ref. |n| by means of the pure-quantum self- consistent harmonic approximation (PQSCHA), 
a semiclassical approach. Data in Fig. |I] refer to g = 5 and different values of damping. The 
reported Monte Carlo data represent the extrapolation to P ^ cxd of the results obtained at 
P = 17, 33, 65, and 129; the reliability of the algorithm is proven by the perfect agreement 
between the exact results and the PIMC data in the non-dissipative system (F = 0). For 
the dissipative model, the PIMC data provide therefore a reference, still lacking until now, 
to verify the validity of the PQSCHA pf. 

Let us turn now to a true many-body dissipative system, and consider the quantum 
discrete 0^ chain, whose Hamiltonian may be written as |l^ 



n 



V{q) 



i=l 



2Rf 



M r 

E 



i=l 



R 



(9) 
(10) 



where v{x) = (1 — x ) /8, Q is the quantum coupling, and and R are the kink energy 
and length, respectively, in the classical continuum limit. 
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FIG. 1: Temperature dependence of the average potential energy {v{x)) for the single particle in 
a quartic double well, for g = 5 and different values of the damping strength T. Empty symbols 
are PIMC data, lines are PQSCHA predictions (Ref. |]ll|) and filled circles are the exact results 
for r = 0. 



Assuming identical independent baths [T^ for each degree of freedom, i.e. Kij{u) 
6ijK{u), Eq. (D gives 

M ^ ^ N 



Zp = Ct^l[ fdqiffl da.kdhk e"^^ 
j=i"^ k=i 



(11) 



where we set qa = + 2 (ciifc cos + bik sin 



2nke 



2itM \ 
P 



again with P = 2N + 1; using the 



dimensionless temperature t = (Pe^) ^, the discretized action reads 



sin — + ^-i^t 



Q'^R P ' 2Rt' 

N 



i=i I k=i 

3R 

+—\l,qi-qi-i)'^+2j2[i(^ik-ai-ik)'^+{hk-bi- 



ik ' ^ik) 



Ik) 



k=l 



+ 



2RtP 



•ii) 



The average quantities for the dissipative 0^ chain presented in Figs. have been 
obtained for periodic chains of length (~ 10^ sites) large enough to be representative of the 
thermodynamic limit for each set of physical parameters and by extrapolating to P ^ cxo 
the results given by simulations at finite P. A Drude-like spectral density was assumed for 
the environmental interaction so that the dissipative kernel reads 

{2TxtklQ) 



° + {2ntk/Q) 



(12) 




FIG. 2: (g?) and {v{qi)) (inset) vs temperature for the (/>^ chain, with Q = 0.2, i? = 5, = 100 
and different values of T. Empty symbols are PIMC data (P ^ oo extrapolations) and lines are 
PQSCHA predictions (Ref. [^). T = 0: circles and solid line; F = 20: squares and short-dashed 
line; F = 100: triangles and long-dashed line. 
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FIG. 3: (qf) vs temperature for the (j)^ chain, with Q = 1, i? = 3, = 10 and different values of 
F. Full symbols are PIMC data at finite Trotter number (P = 81, for finite F, P = 11 for F — > oo); 
empty symbols are P — > oo extrapolations. Lines are PQSCHA predictions. 



where the dissipation strength F and the cut-off frequency are measured in units of 

= QEj^. Comparison of the PIMC results with those of the PQSCHA [|12|, shown in Fig. 
clearly indicates that the predictions of the latter are impressively accurate. Remarkably, 
the accuracy is also preserved for fairly large values of quanticity, quite close to the predicted 
limits of applicability of the PQSCHA scheme, as it appears in Fig. |^. 

The role played by the correction scheme of Eq. (|^) is clearly shown in Fig. ^ (a), where 
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FIG. 4: Internal energy (per site) U for the 4>^ chain, with Q = 1, R = 3, fi^ = 10, and 
different values of T. (a): P-dependence of U for T = 20, and t = 0.2. Full symbols: bare 
PIMC results U{P); empty symbols: harmonically-corrected data C/ha(-P)- Lines are quadratic 
fits, (b): temperature dependence of U. Full symbols: PIMC data for P = 101; empty symbols: 
harmonically-corrected P ^ 00 extrapolations. Lines are PQSCHA predictions. 

^ha(-P) displays a very weak dependence on P, thus allowing to get reliable estimates of 
the P — s> 00 extrapolation even by starting from results obtained for relatively small values 
of P. The relevance of the harmonic correction is more and more evident for increasing 
dissipation strength [Fig. ^ (b)]: indeed, in this regime the finite-P bias essentially comes 
from the bilinear contribution of the influence functional and the latter can be easily healed 
by the above correction scheme, so that, for instance, the characteristic almost flat behavior 
of U at low temperature [r2| is obtained without using exceedingly high values of P. 

In conclusion, we have introduced a formulation of the Path Integral Monte Carlo method 
which greatly simplifies the numerical investigation of the thermodynamics of quantum 
dissipative systems with bilinear influence action, making it possible to afford also many- 
body systems. The combined use of both the discrete Fourier transformed representation 
of the dynamical variables and the correction scheme of Eq. (|^) is essential to efficiently 
circumvent the numerical difficulties arising from the non locality in time and space of the 
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action describing an open many-body system. The application to a double- well potential and 
to a chain have shown the power of the method and have given, as a byproduct, a direct 
confirmation of the validity of the semi-classical approach of Refs. |Tl|, |1^] in the expected 
parameter region, i.e., weak quantum coupling and/or strong dissipation. The numerical 
approach introduced here does not suffer of such limitations, and we expect it to be useful 
for getting meaningful insight on the behavior of strongly quantum systems in presence of 
dissipation, thus making it possible to address interesting problems of mesoscopic physics 
as, e.g., the dissipative transition in Josephson junction arrays |Q. 
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